home *** CD-ROM | disk | FTP | other *** search
/ Pascal Super Library / Pascal Super Library (CW International)(1997).bin / MATH / STAT / STAT.PAS < prev    next >
Pascal/Delphi Source File  |  1990-08-08  |  45KB  |  1,590 lines
  1. {--------------------------------------------------------------------------}
  2. {                         Norton Statistical Library                       }
  3. {                                                                          }
  4. {                              Version   1.00                              }
  5. {                                                                          }
  6. {                                                                          }
  7. {                    Copyright 1990 Norton Associcates                     }
  8. {                           All Rights Reserved                            }
  9. {                          Restricted by License                           }
  10. {--------------------------------------------------------------------------}
  11.  
  12.                      {--------------------------------}
  13.                      {       Unit:   Stat             }
  14.                      {--------------------------------}
  15.  
  16.  
  17. {$S-,R-,V-,D-,A+,B+,N+,E-,I-}
  18.  
  19. UNIT
  20.     stat;
  21.  
  22. INTERFACE
  23.  
  24. CONST
  25.      error_value = -99.9;                         { if any errors }
  26.      maxsize     = 65520;                         { max segment size }
  27.      maxsingle   = maxsize DIV (SIZEOF(SINGLE));  { max element size single array }
  28.      maxdouble   = maxsize DIV (SIZEOF(DOUBLE));  { max element size double array }
  29.      maxlongint  = maxsize DIV (SIZEOF(LONGINT)); { max element size longint array }
  30.  
  31.      v1   : INTEGER = 1;                       { constants for uniform1 }
  32.      v2   : INTEGER = 1000;                    { constants for uniform1 }
  33.      v3   : INTEGER = 30000;                   { constants for uniform1 }
  34.  
  35.      maxorder  = 10;
  36.      maxmatrix = maxorder * 2 - 1;
  37.  
  38.      uno   = 1;
  39.      dos   = 2;
  40.      zero  = 0.0;
  41.      one   = 1.0;
  42.      two   = 2.0;
  43.      c2    = 2.0;
  44.      c3    = -3.0;
  45.      c4    = 4.0;
  46.      c5    = 5.0;
  47.      c6    = 6.0;
  48.      c11   = 11.0;
  49.      c12   = 12.0;
  50.      c17   = 17.0;
  51.      i_4   = 1.0/c4;
  52.      i_16  = 1.0/16.0;
  53.      i_35  = 1.0/35.0;
  54.  
  55. TYPE
  56.  
  57.     single_array_type    = SINGLE;
  58.     single_array_dummy   = ARRAY[1..maxsingle] OF single_array_type;
  59.     single_array_pointer = ^single_array_dummy;
  60.  
  61.     double_array_dummy   = ARRAY[1..maxdouble] OF DOUBLE;
  62.     double_array_pointer = ^double_array_dummy;
  63.  
  64.     longint_array_dummy   = ARRAY[1..maxlongint] OF LONGINT;
  65.     longint_array_pointer = ^longint_array_dummy;
  66.  
  67.     quartype  = ARRAY[1..5] OF SINGLE;
  68.     arry_type = ARRAY[1..maxorder,1..maxorder] OF EXTENDED;
  69.  
  70. PROCEDURE create_single_array( num: WORD; VAR xx : single_array_pointer);
  71. PROCEDURE delete_single_array( num: WORD; VAR xx : single_array_pointer);
  72. PROCEDURE create_longint_array( num: WORD; VAR xx : longint_array_pointer);
  73. PROCEDURE delete_longint_array( num: WORD; VAR xx : longint_array_pointer);
  74.  
  75. FUNCTION uniform1 : SINGLE;
  76. FUNCTION rndnorm1( mean,standev:EXTENDED) :SINGLE;
  77. FUNCTION rndnorm2( mean,standev:EXTENDED) :SINGLE;
  78.  
  79. PROCEDURE insert( n : WORD ; VAR a : single_array_pointer) ;
  80. PROCEDURE qsort( n : WORD; VAR a : single_array_pointer) ;
  81. PROCEDURE remove_avg( n : WORD; VAR a : single_array_pointer ; avg : SINGLE);
  82.  
  83. PROCEDURE means( n :WORD; a : single_array_pointer; VAR xmean,gmean,hmean,rmsmean :SINGLE);
  84. PROCEDURE wxmean( num : WORD; a : single_array_pointer; freq : longint_array_pointer; VAR mean,sd,small,large : SINGLE);
  85. PROCEDURE elem_stat( num:WORD; VAR a:single_array_pointer;
  86.                     VAR small,large:SINGLE; VAR mean,sd:SINGLE);
  87. PROCEDURE moments( n      : WORD;
  88.                    a      : single_array_pointer;
  89.                 VAR ave   : SINGLE;
  90.                 VAR std   : SINGLE;
  91.                 VAR skew  : SINGLE;
  92.                 VAR beta2 : SINGLE);
  93.  
  94. PROCEDURE quart( n : WORD; a : single_array_pointer; VAR quart  : quartype);
  95. FUNCTION percentile( n : WORD; a : single_array_pointer ;percent : SINGLE) : SINGLE;
  96. FUNCTION standard_error(num : WORD ; sd : SINGLE; ntype:WORD) : SINGLE;
  97.  
  98. FUNCTION cdf_prob_to_sd(prob :SINGLE) : SINGLE;
  99. FUNCTION cdf_sd_to_prob(sd:SINGLE) : SINGLE;
  100. FUNCTION int_prob_to_sd(prob :SINGLE) : SINGLE;
  101. FUNCTION int_sd_to_prob(sd:SINGLE) : SINGLE;
  102.  
  103. PROCEDURE corcoef(n: WORD; x,y:single_array_pointer; VAR r:SINGLE);
  104. PROCEDURE autocor(n:WORD; x:single_array_pointer; lag :WORD; VAR auto:single_array_pointer);
  105. PROCEDURE linfit(npts: WORD; x,y,sigmay: single_array_pointer; mode:WORD;VAR a,  b, r:SINGLE);
  106. FUNCTION determ(arry :arry_type; norder : INTEGER) : EXTENDED;
  107. PROCEDURE polfit(npts: WORD;
  108.                 x,y,sdy      : single_array_pointer;
  109.                 nterms,mode  : INTEGER;
  110.                 VAR a        : single_array_pointer;
  111.                 VAR r        : SINGLE;
  112.                 VAR se       : SINGLE);
  113. PROCEDURE mulreg(n : WORD; y,x,z : single_array_pointer;
  114.                           VAR a  : single_array_pointer;
  115.                           VAR r  : SINGLE;
  116.                           VAR se : SINGLE);
  117.  
  118. PROCEDURE smooth121(n:WORD; VAR y: single_array_pointer);
  119. PROCEDURE smooth14641(n:WORD; VAR y: single_array_pointer);
  120. PROCEDURE smoothcurve(n:WORD; VAR y: single_array_pointer);
  121. FUNCTION movavg(n: WORD; a : single_array_pointer; ma:WORD ; k : WORD) : SINGLE;
  122.  
  123. {*****************************************************************************}
  124. {*****************************************************************************}
  125. IMPLEMENTATION
  126. {*****************************************************************************}
  127. {*****************************************************************************}
  128.  
  129. PROCEDURE create_single_array( num: WORD; VAR xx : single_array_pointer);
  130. { Author : Norton Associates
  131.   Purpose: Properly create a single dimensioned heap array
  132.   Version: 1.0
  133.   Date   : 5 May 1990 }
  134.  
  135. VAR
  136.     maxnumber : LONGINT;    { proposed size of array in bytes }
  137.  
  138. BEGIN
  139.      maxnumber := LONGINT(LONGINT(num) * LONGINT(SIZEOF(single_array_type)));
  140.      IF (num < uno) or (num > maxsingle) THEN
  141.      BEGIN
  142.           WRITELN('Sorry the single array size is to large to create ');
  143.           WRITELN('You wanted = ',num:10,', The max single size = ',maxsingle:10);
  144.           HALT;
  145.      END;
  146.      GETMEM(xx,maxnumber);
  147. END;
  148.  
  149. PROCEDURE delete_single_array( num: WORD; VAR xx : single_array_pointer);
  150. { Author : Norton Associates
  151.   Purpose: Properly delete a single dimensioned heap array
  152.   Version: 1.0
  153.   Date   : 5 May 1990 }
  154.  
  155. VAR
  156.    maxnumber : LONGINT;  { proposed size of array in bytes }
  157.  
  158. BEGIN
  159.      maxnumber := LONGINT(LONGINT(num) * LONGINT(SIZEOF(single_array_type)));
  160.      IF maxnumber > maxsize THEN
  161.      BEGIN
  162.           WRITELN('sorry the single array size is to large to delete ',maxnumber);
  163.           HALT;
  164.      END;
  165.      FREEMEM(xx,maxnumber);
  166. END;
  167.  
  168. PROCEDURE create_longint_array( num: WORD; VAR xx : longint_array_pointer);
  169. { Author : Norton Associates
  170.   Purpose: Properly create a longint dimensioned heap array
  171.   Version: 1.0
  172.   Date   : 5 May 1990 }
  173.  
  174. VAR
  175.    maxnumber : LONGINT;  { proposed size of array in bytes }
  176.  
  177. BEGIN
  178.      maxnumber := LONGINT(LONGINT(num) * LONGINT(SIZEOF(single_array_type)));
  179.      IF (num < uno) or (num > maxlongint) THEN
  180.      BEGIN
  181.           WRITELN('sorry longint array size is to large to create',maxlongint);
  182.           HALT;
  183.      END;
  184.      GETMEM(xx,maxnumber);
  185. END;
  186.  
  187. PROCEDURE delete_longint_array( num: WORD; VAR xx : longint_array_pointer);
  188. { Author : Norton Associates
  189.   Purpose: Properly delete a longint dimensioned heap array
  190.   Version: 1.0
  191.   Date   : 5 May 1990 }
  192.  
  193. VAR
  194.    maxnumber : LONGINT;  { proposed size of array in bytes }
  195.  
  196. BEGIN
  197.      maxnumber := LONGINT(LONGINT(num) * LONGINT(SIZEOF(single_array_type)));
  198.      IF maxnumber > maxsize THEN
  199.      BEGIN
  200.           WRITELN('sorry longint array size is to large to delete',maxnumber);
  201.           HALT;
  202.      END;
  203.      GETMEM(xx,maxnumber);
  204. END;
  205.  
  206. PROCEDURE remove_avg(n : WORD; VAR a : single_array_pointer ; avg : SINGLE);
  207. { Author : Norton Associates
  208.   Purpose: Remove a constant value from an array
  209.   Version: 1.0
  210.   Date   : 5 May 1990 }
  211.  
  212. VAR
  213.    j : WORD;
  214.  
  215. BEGIN
  216.    IF n > 0 THEN
  217.    BEGIN
  218.        FOR j := uno TO n DO
  219.           a^[j] := a^[j] - avg;
  220.    END
  221.    ELSE
  222.        WRITELN('remavg> n < 1: No data');
  223. END;
  224.  
  225. PROCEDURE wxmean(num : WORD; a : single_array_pointer; freq : longint_array_pointer; VAR mean,sd,small,large : SINGLE);
  226. { Author : Norton Associates
  227.   Purpose: Calculate a weighted arithemetic mean
  228.   Version: 1.0
  229.   Date   : 5 May 1990 }
  230.  
  231. VAR
  232.     number,                 { num = number of variables to process }
  233.     j        : WORD;        { loop index }
  234.     dum,
  235.     sum,sum2 : EXTENDED;    { sum of squares }
  236.  
  237. BEGIN
  238.  
  239. { check to see if there is enough data }
  240.     IF num < dos THEN
  241.     BEGIN
  242.         WRITELN('Wxmean> No sample data: n < 2');
  243.         small := error_value;
  244.         large := error_value;
  245.         mean  := error_value;
  246.         sd    := error_value;
  247.         EXIT;
  248.     END;
  249.  
  250. { zero out and initalize some data }
  251.     sum   := zero;
  252.     sum2  := zero;
  253.     small := a^[1];
  254.     large := a^[1];
  255.  
  256. { calculate sum of squares }
  257.     FOR j := uno TO num DO
  258.     BEGIN
  259.         dum  := a^[j] * freq^[j];
  260.         sum  := sum + dum;
  261.         sum2 := sum2 + SQR(dum);
  262.         number := number + freq^[j];
  263.         IF a^[j] > large THEN
  264.            large := a^[j]
  265.         ELSE IF a^[j] < small THEN
  266.                 small := a^[j];
  267.     END;
  268.  
  269. { calculate mean of sample population }
  270.     mean := sum / number;
  271.  
  272. { check to see if standard deviation can be calculated }
  273.     IF number - uno < uno THEN
  274.     BEGIN
  275.        sd := error_value;
  276.        WRITELN('wxmean> no standard deviation calculated');
  277.        EXIT;
  278.     END
  279.     ELSE IF (sum2 - (number * mean * mean) ) <= zero THEN
  280.             sd := zero
  281.          ELSE
  282.             sd := SQRT((sum2 - (number * mean * mean) ) / (number-1));
  283. END;
  284.  
  285. PROCEDURE means(n :WORD; a : single_array_pointer; VAR xmean,gmean,hmean,rmsmean:SINGLE);
  286. { Author : Norton Associates
  287.   Purpose: Calculate the arithemetic, geometric, root mean square and
  288.            harmonic mean
  289.   Version: 1.0
  290.   Date   : 5 May 1990 }
  291.  
  292. VAR
  293.    xsum,gsum,hsum,rmssum : EXTENDED;
  294.    gboolean,hboolean     : BOOLEAN;
  295.    j                     : WORD;
  296. BEGIN
  297.  
  298. { clear and set up the data }
  299.      xsum := zero;                { arithmetic mean }
  300.      hsum := zero;                { harmonic mean }
  301.      gsum := zero;                { geometric mean }
  302.      rmssum   := zero;            { roor mean square mena }
  303.      gboolean := TRUE;
  304.      hboolean := TRUE;
  305.  
  306. { check to see if there is enough data }
  307.      IF n  < dos THEN
  308.      BEGIN
  309.           WRITELN('means> no sample data: n < 2');
  310.           EXIT;
  311.      END;
  312.  
  313. { sort out the data }
  314.      FOR j := uno TO n DO
  315.      BEGIN
  316.          xsum := xsum + a^[j];
  317.          IF a^[j] > zero THEN
  318.             gsum := gsum + LN(a^[j])
  319.          ELSE
  320.             gboolean := FALSE;              { bad data flag }
  321.          IF a^[j] > zero THEN
  322.             hsum := hsum + one / a^[j]
  323.          ELSE
  324.             hboolean := FALSE;              { bad data flag }
  325.          rmssum := rmssum + SQR(a^[j]);
  326.      END;
  327.  
  328.      IF gboolean = FALSE THEN
  329.      BEGIN
  330.         WRITELN('means> Gmean can not be calculated: Hmean = -99.9');
  331.         gmean := error_value
  332.      END
  333.      ELSE
  334.         gmean := EXP(gsum/n);
  335.      IF hboolean = FALSE THEN
  336.      BEGIN
  337.         WRITELN('means> Hmean can not be calculated: Hmean = -99.9');
  338.         hmean := error_value
  339.      END
  340.      ELSE
  341.         gmean := EXP(gsum/n);
  342.      xmean := xsum / n;
  343.      rmsmean := SQRT(rmssum / n);
  344. END;
  345.  
  346. FUNCTION cdf_prob_to_sd(prob :SINGLE) : SINGLE;
  347. { Author : Norton Associates
  348.   Purpose: Calculate the standard deviation from probability
  349.   Version: 1.0
  350.   Date   : 5 May 1990
  351.   PROB=PROBABILITY (0-1) INPUT
  352.   SD=NO. OF STANDARD DEVIATIONS,OUTPUT  }
  353.  
  354. CONST
  355.         c0 : extended = 2.515517;
  356.         c1 : extended = 0.802853;
  357.         c2 : extended = 0.010328;
  358.         d1 : extended = 1.432788;
  359.         d2 : extended = 0.189269;
  360.         d3 : extended = 0.001308;
  361. VAR
  362.         dn,up,xp : EXTENDED;
  363.         t1,t2,t3 : EXTENDED;
  364.  
  365. BEGIN
  366.  
  367. { check to see if the probability ( prob ) is in range }
  368.      IF (prob < zero) OR (prob > one) THEN
  369.      BEGIN
  370.          WRITELN('prob_to_sd> input probability out of range (0-1)');
  371.          EXIT;
  372.      END;
  373.      IF prob > 0.5 THEN
  374.           xp := one - prob
  375.      ELSE
  376.           xp := prob;
  377.      t1 := SQRT(LN(one/SQR(xp)));
  378.      t2 := t1 * t1;
  379.      t3 := t2 * t1;
  380.      up := c0  + c1*t1 + c2*t2;
  381.      dn := one + d1*t1 + d2*t2 + d3*t3;
  382.      xp := t1 - (up/dn);
  383.      IF prob <= 0.5 THEN
  384.         cdf_prob_to_sd := xp * (-one);
  385.      cdf_prob_to_sd := xp;
  386. END;
  387.  
  388. FUNCTION cdf_sd_to_prob(sd:SINGLE) : SINGLE;
  389. { Author : Norton Associates
  390.   Purpose: Calculate the standard deviation from probability
  391.   Version: 1.0
  392.   Date   : 5 May 1990
  393.   sd = (mean-a)/standard deviation
  394.   prob = percentile of a between 0 and 100%)
  395. }
  396.  
  397. CONST
  398.         p  : extended =  0.231641900;
  399.         b1 : extended =  0.319381530;
  400.         b2 : extended = -0.356563782;
  401.         b3 : extended =  1.781477937;
  402.         b4 : extended = -1.821255978;
  403.         b5 : extended =  1.330274429;
  404. VAR
  405.    x,y1,y2,y3,y4,y5    : EXTENDED;
  406.    t,z,r               : EXTENDED;
  407.  
  408. BEGIN
  409.  
  410. { check to see if the standard deviation ( sd ) is in range }
  411.       IF sd < zero THEN
  412.       BEGIN
  413.          WRITELN('cdf_sd_prob> input standard deviation out of range: sd >=0.0');
  414.          EXIT;
  415.       END;
  416.       x := sd;
  417.       r := EXP(-(x*x)/two)/2.5066282746;
  418. {  r = frequency }
  419.       z  := x;
  420.       y1 := one/(one+(p*ABS(x)));
  421.       y2 := y1 * y1;
  422.       y3 := y2 * y1;
  423.       y4 := y3 * y1;
  424.       y5 := y4 * y1;
  425.       t  := one - r*(b1*y1+ b2*y2 + b3*y3 + b4*y4 + b5*y5);
  426.       IF z > zero THEN
  427.          cdf_sd_to_prob := t
  428.       ELSE
  429.          cdf_sd_to_prob := one - t;
  430. { prob = [0% to 100%]}
  431. END;
  432.  
  433. FUNCTION int_prob_to_sd(prob :SINGLE) : SINGLE;
  434. { Author : Norton Associates
  435.   Purpose: Calculate the probability from the standard deviation
  436.   Version: 1.0
  437.   Date   : 5 May 1990 }
  438.  
  439. VAR
  440.    dummy : SINGLE;
  441.  
  442. BEGIN
  443.  
  444. { check to see if the probability is in range }
  445.      IF (prob < zero) OR (prob > one) THEN
  446.      BEGIN
  447.         WRITELN('int_prob_to_sd> input probability out of range: prob (0-1)');
  448.         EXIT;
  449.      END;
  450.      dummy := 0.50 + (0.50 * prob) ;
  451.      int_prob_to_sd := cdf_prob_to_sd(dummy);
  452. END;
  453.  
  454. FUNCTION int_sd_to_prob(sd:SINGLE) : SINGLE;
  455. { Author : Norton Associates
  456.   Purpose: Calculate the probability from the standard deviation
  457.   Version: 1.0
  458.   Date   : 5 May 1990 }
  459.  
  460. VAR
  461.    dummy: SINGLE;
  462.  
  463. BEGIN
  464.  
  465. { check to see if the standard deviation ( sd ) is in range }
  466.   IF sd <= zero THEN
  467.   BEGIN
  468.     WRITELN('int_sd_to_prob> negative or zero standard deviation not allowed');
  469.     EXIT;
  470.   END;
  471.   dummy := cdf_sd_to_prob(sd);
  472.   int_sd_to_prob := dummy - (0.50 - (dummy - 0.50) )
  473. END;
  474.  
  475. PROCEDURE linfit(npts : WORD; x,y,sigmay: single_array_pointer;
  476.                       mode : WORD; VAR a,  b, r : SINGLE);
  477. { Author : Norton Associates
  478.   Purpose: Calculate a linear least squares fit for x and y pairs of data
  479.   Version: 1.0
  480.   Date   : 5 May 1990 }
  481.  
  482. VAR
  483.   sum,sumx,sumy,sumx2,sumy2,sumxy : EXTENDED;
  484.   x1,y1,weight,delta              : EXTENDED;
  485.   i                               : WORD;
  486. { x,y pairs of data
  487.   sigmay standard deviation of y
  488.   npts  =  number of pairs of as
  489.   mode = 0 =  no weighting
  490.   mode = 1 =  instrumental weighting (sigmay array must exist)
  491.      a  =  intercept
  492.      b  =  slope
  493.      r  =  correlation coefficient }
  494.  
  495. BEGIN
  496.  
  497. { zero out items for use }
  498.       sum    :=  zero;
  499.       sumx   :=  zero;
  500.       sumy   :=  zero;
  501.       sumxy  :=  zero;
  502.       sumx2  :=  zero;
  503.       sumy2  :=  zero;
  504. { sum over npts pairs of points }
  505.  
  506. { check to see if there is enough data }
  507.       IF npts < dos THEN
  508.       BEGIN
  509.            WRITELN('linfit> number of data points < 2');
  510.            EXIT;
  511.       END;
  512.  
  513. { calculate sum of square for the data }
  514.       FOR i  :=  uno TO npts DO
  515.       BEGIN
  516.          x1  :=  x^[i];
  517.          y1  :=  y^[i];
  518.          IF (mode  = 0) OR (y1 = zero) THEN
  519.             weight := one
  520.          ELSE
  521.             weight  :=  one / SQR(sigmay^[i]);
  522.          sum    :=  sum   + weight;
  523.          sumx   :=  sumx  + weight * x1;
  524.          sumy   :=  sumy  + weight * y1;
  525.          sumx2  :=  sumx2 + weight * x1 * x1;
  526.          sumy2  :=  sumy2 + weight * y1 * y1;
  527.          sumxy  :=  sumxy + weight * x1 * y1;
  528.       END;
  529. { calculate final answers }
  530.  
  531.    delta  :=  sum * sumx2 - sumx * sumx;
  532.    a      :=  (sumx2*sumy - sumx*sumxy) / delta;
  533.    b      :=  (sumxy*sum  - sumx*sumy) / delta;
  534.    r      :=  (sum*sumxy  - sumx*sumy) /
  535.               SQRT(delta  * (sum * sumy2 - sumy * sumy));
  536. END;
  537.  
  538. FUNCTION standard_error(num : WORD ; sd : SINGLE; ntype : WORD) : SINGLE;
  539. { Author : Norton Associates
  540.   Purpose: Calculate the standard error based upon a normal distribution
  541.   Version: 1.0
  542.   Date   : 5 May 1990 }
  543.  
  544. BEGIN
  545.  
  546. { check to see if there is enough data }
  547.   IF num < dos THEN
  548.   BEGIN
  549.        WRITELN('Standard_error> number of samples < 2');
  550.        EXIT;
  551.   END;
  552.  
  553. { determine the type and calculate the standard error }
  554.   CASE ntype OF
  555.      1 : standard_error := sd / SQRT(num);
  556.      2 : standard_error := sd * SQRT(num);
  557.      3 : standard_error := 1.2533 * sd / SQRT(num);
  558.      4 : standard_error := sd / SQRT(two * num);
  559.      5 : standard_error := 0.6028 * sd / SQRT(num);
  560.      6 : standard_error := sd * SQRT(one + (two * SQR(sd))) / SQRT(two * num);
  561.      7 : standard_error := (one - SQR(sd)) / SQRT(num);
  562.   END;
  563. END;
  564.  
  565. PROCEDURE qsort(    n   : WORD     ;   { number of as }
  566.                 VAR a   : single_array_pointer) ;   { number of as in a }
  567. { Author : Norton Associates
  568.   Purpose: Sort data based upon quick sort algorithm
  569.   Version: 1.0
  570.   Date   : 5 May 1990 }
  571.  
  572. CONST
  573.      brkeven  = 13; { more than 10 as use qsort otherwise use insert sort}
  574.      maxstack = 20; { enough stack space for 1_000_000 random numbers }
  575.  
  576. VAR
  577.      low,high,ti,median,jd,
  578.      bi,id,index,pt,j,i    : WORD;
  579.      stack                 : 0..maxstack;
  580.      top,bottom            : ARRAY[1..maxstack] OF WORD;
  581.      pivot                 : SINGLE;   { change to match singlearray }
  582.  
  583. PROCEDURE swap_position( median,low : WORD);
  584.  
  585. VAR
  586.    t : SINGLE;
  587.  
  588. BEGIN
  589.     t          := a^[median];
  590.     a^[median] := a^[low];
  591.     a^[low]    := t;
  592. END;
  593.  
  594. BEGIN { main program }
  595.  
  596.      { non-recursive qsort, must use stack : 10% performance increase }
  597.  
  598.      { initialize the qsort stack }
  599.      stack         := uno;
  600.      top[stack]    := uno;
  601.      bottom[stack] := n;
  602.  
  603. { check to see if there is enough data }
  604.      IF n < dos THEN
  605.      BEGIN
  606.           WRITELN('qsort> no sample data');
  607.           EXIT;
  608.      END;
  609.      REPEAT
  610.        { pop information off the stack }
  611.            low  := top[stack];
  612.            high := bottom[stack];
  613.  
  614.      { use of insert sort for small sublists : 16% performance increase }
  615.      { use qsort for more than brkeven or use insert sort for small sublists }
  616.            WHILE (high - low) > brkeven DO
  617.            BEGIN
  618.                ti := low;
  619.                bi := high;
  620.    { median of three rule : 10-100% performance increase, almost no worst case }
  621.                median := (low + high) DIV 2;   { initial guess }
  622.                IF a^[low] > a^[median] THEN
  623.                     swap_position(median,low);
  624.  
  625.                IF a^[high] < a^[median]  THEN
  626.                BEGIN
  627.                     swap_position(high,median);
  628.                     IF a^[low] > a^[median] THEN
  629.                          swap_position(low,median);
  630.                END;
  631.  
  632.                pivot := a^[median];  { pivot = median of three }
  633.  
  634.         { the old qsort partition algorithm, median of three optimized}
  635.                REPEAT
  636.                      REPEAT dec(bi); UNTIL a^[bi] <= pivot;
  637.                      REPEAT inc(ti); UNTIL a^[ti] >= pivot;
  638.                      IF ti <= bi THEN
  639.                         swap_position(bi,ti);
  640.                UNTIL ti > bi ;
  641.  
  642.         { larger of the two stacks saved first: 10% performance increase}
  643.                IF (bi - low) > (high - ti) THEN
  644.                BEGIN
  645.                    top[stack]    := low;
  646.                    bottom[stack] := bi;
  647.                    low           := ti;
  648.                END
  649.                ELSE
  650.                BEGIN
  651.                    top[stack]    := ti;
  652.                    bottom[stack] := high;
  653.                    high          := bi;
  654.                END;
  655.         { wrap up non-recursive part of qsort }
  656.                inc(stack);
  657.                IF stack > maxstack  THEN
  658.                BEGIN
  659.                     WRITELN('qsort> stack space exceeded');
  660.                     HALT;
  661.                END;
  662.            END;
  663.  
  664.      { use insertion sort for small sublists : 16% performance increase}
  665.            IF (high - low) > 0 THEN
  666.            BEGIN
  667.              id := low + uno;
  668.              FOR i := id TO high DO
  669.              BEGIN
  670.                   pivot := a^[i];
  671.                   IF a^[low] > pivot THEN
  672.                   BEGIN
  673.                        FOR j := i DOWNTO id DO
  674.                            a^[j] := a^[j-1];
  675.                        a^[low] := pivot;
  676.                   END
  677.                   ELSE
  678.                   BEGIN
  679.                        j := i;
  680.                        WHILE a^[j-1] > pivot DO
  681.                        BEGIN
  682.                            a^[j] := a^[j-1];
  683.                            dec(j);
  684.                        END;
  685.                        a^[j] := pivot;
  686.                   END
  687.              END;
  688.            END;
  689.            dec(stack);
  690.      UNTIL stack = 0
  691. END;
  692.  
  693. PROCEDURE insert(n   : WORD     ;   { number of as }
  694.                 VAR a   : single_array_pointer) ;   { number of as in a }
  695. { Author : Norton Associates
  696.   Purpose: Sort data based upon insertion sort algorithm
  697.   Version: 1.0
  698.   Date   : 5 May 1990 }
  699.  
  700. VAR
  701.    i,j   : WORD;
  702.    pivot : SINGLE;
  703.  
  704. BEGIN
  705.  
  706. { check to see if there is enough data }
  707.      IF n < dos THEN
  708.      BEGIN
  709.           WRITELN('insert> no sample data');
  710.           EXIT;
  711.      END;
  712.      FOR i := 2 TO n DO
  713.      BEGIN
  714.           pivot := a^[i];
  715.           IF a^[1] > pivot THEN
  716.           BEGIN
  717.                FOR j := i DOWNTO 2 DO
  718.                    a^[j] := a^[j-1];
  719.                a^[1] := pivot;
  720.           END
  721.           ELSE
  722.           BEGIN
  723.                j := i;
  724.                WHILE a^[j-1] > pivot DO
  725.                BEGIN
  726.                     a^[j] := a^[j-1];
  727.                     dec(j);
  728.                END;
  729.                a^[j] := pivot;
  730.           END
  731.      END;
  732. END;
  733.  
  734. FUNCTION uniform1 : SINGLE;
  735. { Author : Norton Associates
  736.   Purpose: Sort data based upon insertion sort algorithm
  737.   Version: 1.0
  738.   Date   : March 1987 Byte Magazine }
  739.  
  740. VAR
  741.    temp : SINGLE ;
  742.  
  743. BEGIN
  744.    v1 := 171  * ( v1 MOD 177) - 2 * ( v1 DIV 177);
  745.    IF v1 < 0 THEN
  746.       v1 := v1 + 30269;
  747.    v2 := 172 * (v2 MOD 176) - 35 * ( v2 DIV 176);
  748.    IF v2 < 0 THEN
  749.       v2 := v2 + 30307;
  750.    v3 := 170 * (v3 MOD 178) - 63 * ( v3 MOD 178);
  751.    IF v3 < 0 THEN
  752.       v3 := v3 + 30323;
  753.    temp := (v1 / 30269.0) + (v2 / 30307.0) + (v3 / 30323.0);
  754.    uniform1 := temp - TRUNC(temp);
  755. END;
  756.  
  757.  
  758. FUNCTION rndnorm1(mean , standev : EXTENDED) :SINGLE;
  759. { Author : Norton Associates
  760.   Purpose: Calculate random normal number
  761.   Version: 1.0
  762.   Date   : 5 May 1990 }
  763.  
  764. VAR
  765.    randoma,randomb,randomc,radius2,deviate  : EXTENDED;
  766.  
  767. BEGIN
  768.      REPEAT
  769.            randoma := two * RANDOM - one;
  770.            randomb := two * RANDOM - one;
  771.            radius2 := SQR(randoma) + SQR(randomb);
  772.      UNTIL radius2 < one;
  773.  
  774.      IF RANDOM(2) = uno THEN
  775.         randomc := randomb
  776.      ELSE
  777.         randomc := randoma;
  778.  
  779.      deviate := randomc * SQRT(( - two * LN(radius2)) / radius2);
  780.      rndnorm1 := mean + deviate * standev;
  781. END;
  782.  
  783. FUNCTION rndnorm2( mean , standev : EXTENDED) :SINGLE;
  784. { Author : Norton Associates
  785.   Purpose: Calculate random normal number ( approximation  to rndnom1 )
  786.   Version: 1.0
  787.   Date   : 5 May 1990 }
  788.  
  789. VAR
  790.    sum : SINGLE;
  791.    j   : WORD;
  792.  
  793. BEGIN
  794.     sum   := zero;
  795.     FOR j := uno TO 12 DO
  796.         sum := sum + RANDOM;
  797.     rndnorm2 := mean + (sum - 6.0) * standev;
  798. END;
  799.  
  800. FUNCTION power( number, exponent : SINGLE) : SINGLE;
  801. { Author : Norton Associates
  802.   Purpose: Raise a number to a number
  803.   Version: 1.0
  804.   Date   : 5 May 1990 }
  805.  
  806. BEGIN
  807.      IF (number <> zero) AND (exponent <> zero) THEN
  808.         power := EXP( exponent * LN(number) )
  809.      ELSE
  810.          WRITELN('power> can not take power of ',number:10,exponent:10);
  811. END;
  812.  
  813.  
  814. PROCEDURE elem_stat(num:WORD; VAR a:single_array_pointer;
  815.                VAR small,large:SINGLE; VAR mean,sd:SINGLE);
  816. { Author : Norton Associates
  817.   Purpose: Calculate mean, standard deviation the largest and smallest value
  818.   Version: 1.0
  819.   Date   : 5 May 1990 }
  820.  
  821. VAR
  822.     j        : WORD;        { loop index }
  823.     sum,sum2 : EXTENDED;    { sum of squares }
  824.  
  825. BEGIN
  826.  
  827. { check to see if there is data }
  828.     IF num < uno THEN
  829.     BEGIN
  830.         WRITELN('elem_stat> no data points');
  831.         small := error_value;
  832.         large := error_value;
  833.         mean  := error_value;
  834.         sd    := error_value;
  835.         EXIT;
  836.     END;
  837.  
  838. { zero out data }
  839.     sum   := zero;
  840.     sum2  := zero;
  841.     small := a^[1];
  842.     large := a^[1];
  843.  
  844. { calculate sum of squares }
  845.     FOR j := uno TO num DO
  846.     BEGIN
  847.         sum  := sum + a^[j];
  848.         sum2 := sum2 + SQR(a^[j]);
  849.         IF a^[j] > large THEN
  850.            large := a^[j]
  851.         ELSE IF a^[j] < small THEN
  852.                 small := a^[j];
  853.     END;
  854.  
  855. { calculate mean and standard deviation }
  856.     mean := sum / num;
  857.     IF num - uno < uno THEN
  858.     BEGIN
  859.        sd := error_value;
  860.        WRITELN('elem_stat> number of samples = 1 no standard deviation');
  861.        EXIT;
  862.     END
  863.     ELSE IF (sum2 - (num * mean * mean) ) <= zero THEN
  864.             sd := zero
  865.          ELSE
  866.             sd := SQRT((sum2 - (num * mean * mean) ) / (num-1));
  867. END;
  868.  
  869. PROCEDURE moments( n      : WORD;
  870.                    a      : single_array_pointer;
  871.                 VAR ave   : SINGLE;
  872.                 VAR std   : SINGLE;
  873.                 VAR skew  : SINGLE;
  874.                 VAR beta2 : SINGLE);
  875. { Author : Norton Associates
  876.   Purpose: Calculate the first 4 moments of the sample population
  877.   Version: 1.0
  878.   Date   : 5 May 1990 }
  879. { COMPUTE STATISTICAL MOMENTS  }
  880.  
  881. VAR
  882.    d,sum,sum2,sum3,sum4,z2,z3,z4,s2,s3,s4 : EXTENDED;
  883.    t : SINGLE;
  884.    i : WORD;
  885.  
  886. BEGIN
  887.  
  888. { check to see if there is enough data }
  889.       IF n < dos THEN
  890.       BEGIN
  891.          WRITELN('moments> no sample data');
  892.          EXIT;
  893.       END;
  894. {...COMPUTE MEAN }
  895.       t   :=  n;
  896.       sum :=  zero;
  897.       FOR i  :=  uno TO n DO
  898.           sum :=  sum + a^[i];
  899.       ave :=  sum / t;
  900. {...COMPUTE OTHER MOMENTS }
  901.       sum2 :=  zero;
  902.       sum3 :=  zero;
  903.       sum4 :=  zero;
  904.       FOR  i := uno TO n DO
  905.       BEGIN
  906.          d  := a^[i] - ave;
  907.          z2 := d * d;
  908.          z3 := d * z2;
  909.          z4 := d * z3;
  910.          sum2 := sum2 + z2;
  911.          sum3 := sum3 + z3;
  912.          sum4 := sum4 + z4;
  913.      END;
  914. {...COMPUTE VARIANCE,STANDARD DEVIATION   }
  915.       s2  := sum2 /(t-one);
  916.       std := SQRT(s2);
  917. {...COMPUTE COEFFICIENT OF SKEWNESS  }
  918.      IF s2 <> zero THEN
  919.       BEGIN
  920.          s3 := sum3 / t;
  921.          skew := s3 / power(s2,1.5);
  922. {...COMPUTE COEFFICIENT OF KURTOSIS   }
  923.          s4 := sum4 / t;
  924.          beta2 := s4 / (s2*s2);
  925.       END
  926.       ELSE
  927.       BEGIN
  928.           WRITELN('moments> no data for beta2 and skewness');
  929.           beta2 := error_value;
  930.           skew  := error_value;
  931.       END;
  932. END;
  933.  
  934. PROCEDURE quart(    n      : WORD;
  935.                     a  : single_array_pointer;
  936.                 VAR quart  : quartype);
  937. { Author : Norton Associates
  938.   Purpose: Calculate the three quartiles and the smallest and largest
  939.            of sample population
  940.   Version: 1.0
  941.   Date   : 5 May 1990 }
  942.  
  943. CONST
  944.      p : ARRAY[1..3] OF SINGLE = (0.25,0.50,0.75);
  945.  
  946. VAR
  947.      r,add,dn,dr :SINGLE;
  948.      nt,nb,j     : WORD;
  949. BEGIN
  950.  
  951. { check to see if there is data }
  952.      IF n < dos THEN
  953.      BEGIN
  954.           WRITELN('quart> no sample data');
  955.           EXIT;
  956.      END;
  957.      qsort(n,a);
  958.      quart[1] := a^[1];         { smallest }
  959.      quart[5] := a^[n];         { largest }
  960.      FOR j := uno TO 3 DO
  961.      BEGIN
  962.         r := p[j] * n;
  963.         nt := ROUND(r + 0.999);
  964.         nb := TRUNC(r);
  965.         add := zero;
  966.         IF (nb = 0) OR (nt = 0) THEN
  967.            quart[j+1] := a^[1]
  968.         ELSE
  969.         BEGIN
  970.            IF nb <> nt THEN
  971.            BEGIN
  972.              dn  := r - nb;
  973.              dr  := a^[nt] - a^[nb];
  974.              add := dn * dr;
  975.            END;
  976.            quart[j+1] := a^[nb] + add;
  977.         END;
  978.      END;
  979. END;
  980.  
  981. FUNCTION percentile(n : WORD; a : single_array_pointer ; percent : SINGLE) : SINGLE;
  982. { Author : Norton Associates
  983.   Purpose: Calculate the value in the sample population based upon
  984.            percentile provided
  985.   Version: 1.0
  986.   Date   : 5 May 1990 }
  987.  
  988. VAR
  989.      r,add,dn,dr :SINGLE;
  990.      nt,nb,j     : WORD;
  991.  
  992. BEGIN
  993.  
  994. { check to see if percentile is in range }
  995.      IF percent < zero  THEN
  996.      BEGIN
  997.          WRITELN('percentile> error percentile must be >= 0.0');
  998.          EXIT;
  999.      END;
  1000.  
  1001. { any data }
  1002.      IF n < one  THEN
  1003.      BEGIN
  1004.          WRITELN('percentile> no sample data');
  1005.          EXIT;
  1006.      END;
  1007.  
  1008.      percent := percent / 100.0;
  1009.      qsort(n,a);
  1010.      r   := percent * n;
  1011.      nt  := ROUND(r + 0.999);
  1012.      nb  := TRUNC(r);
  1013.      add := zero;
  1014.      IF (nb = 0) OR (nt = 0) THEN
  1015.         percentile := a^[1]
  1016.      ELSE
  1017.      BEGIN
  1018.         IF nb <> nt THEN
  1019.         BEGIN
  1020.            dn  := r - nb;
  1021.            dr  := a^[nt] - a^[nb];
  1022.            add := dn * dr;
  1023.         END;
  1024.         percentile := a^[nb] + add;
  1025.      END;
  1026. END;
  1027.  
  1028. FUNCTION determ(arry :arry_type; norder : INTEGER) : EXTENDED;
  1029. { Author : Norton Associates
  1030.   Purpose: Solve simultaneous equations
  1031.   Version: 1.0
  1032.   Date   : 5 May 1990 }
  1033.  
  1034. VAR
  1035.       det,save : EXTENDED;
  1036.       j,k,k1,
  1037.       l,i      : WORD;
  1038.  
  1039. BEGIN
  1040.    det := one;
  1041.  
  1042.    FOR k := uno TO norder DO
  1043.    BEGIN
  1044.       IF arry[k,k] = zero THEN
  1045.       BEGIN
  1046.          FOR j := k TO norder DO
  1047.          BEGIN
  1048.             IF arry[k,j] = zero  THEN
  1049.             BEGIN
  1050.                  WRITELN('determinate> determinate contains a zero');
  1051.                  determ := zero;
  1052.                  EXIT;
  1053.             END;
  1054.          END;
  1055.  
  1056.          FOR i := k TO norder DO
  1057.          BEGIN
  1058.             save      := arry[i,j];
  1059.             arry[i,j] := arry[i,k];
  1060.             arry[i,k] := save;
  1061.          END;
  1062.          det := - det;
  1063.       END;
  1064.       det := det * arry[k,k];
  1065.       IF k - norder < 0 THEN
  1066.       BEGIN
  1067.          k1 := k + uno;
  1068.          FOR i := k1 TO norder DO
  1069.             FOR j := k1 TO norder DO
  1070.                arry[i,j] := arry[i,j] - arry[i,k] * arry[k,j] / arry[k,k];
  1071.       END;
  1072.    END;
  1073.    determ := det;
  1074. END;
  1075.  
  1076. PROCEDURE polfit(npts       : WORD;
  1077.                 x,y,sdy     : single_array_pointer;
  1078.                 nterms,mode : INTEGER;
  1079.                 VAR a       : single_array_pointer;
  1080.                 VAR r       : SINGLE;
  1081.                 VAR se      : SINGLE);
  1082. { Author : Norton Associates
  1083.   Purpose: Determine from least squares the regression
  1084.            and correlation coefficients
  1085.   Version: 1.0
  1086.   Date   : 5 May 1990 }
  1087.  
  1088. VAR
  1089.   wt,xterm,yterm,
  1090.   delta,srs,sum_y,sum_y2,
  1091.   ycalc,x_power,residual : EXTENDED;
  1092.   sumx                   : ARRAY[1..maxmatrix] OF EXTENDED;
  1093.   sumy                   : ARRAY[1..maxorder] OF EXTENDED;
  1094.   arry                   : arry_type;
  1095.   n,j,nmax,k,l           : WORD;
  1096.   xj,yj                  : SINGLE;
  1097.   sigmay2,se2,dummy      : EXTENDED;
  1098.  
  1099. PROCEDURE bad_data(nn : WORD);
  1100. VAR
  1101.    j : WORD;
  1102. BEGIN
  1103.      r      := zero;
  1104.      se     := zero;
  1105.      FOR j := uno TO nterms DO
  1106.          a^[j] := zero;
  1107.      WRITELN('polfit> Determinate or correlation approximately 0.0   ',nn);
  1108. END;
  1109.  
  1110. BEGIN
  1111.  
  1112. { check to see if there is data }
  1113.       IF npts < dos THEN
  1114.       BEGIN
  1115.            WRITELN('polfit> number of data points < 2');
  1116.            EXIT;
  1117.       END;
  1118.  
  1119. { check the number of terms versus maximum order }
  1120.       IF nterms > maxorder THEN
  1121.       BEGIN
  1122.          WRITELN('polfit> maxorder for polfit is ',maxorder:10);
  1123.          EXIT;
  1124.       END;
  1125.  
  1126. { check to see if number of points is more than number of terms }
  1127.       IF nterms >= npts THEN
  1128.       BEGIN
  1129.          WRITELN('polfit> number of terms must be less than number of points');
  1130.          EXIT;
  1131.       END;
  1132. { zero out data }
  1133.       r := zero;
  1134.       FOR j := uno TO nterms DO
  1135.           a^[j] := zero;
  1136.       nmax := 2 * nterms - uno;
  1137.       FOR n := uno TO nmax DO
  1138.           sumx[n] := zero;
  1139.  
  1140. { calculate sum of squares for problem }
  1141.       FOR j := uno TO nterms DO
  1142.           sumy[j] := zero;
  1143.       FOR j := uno TO npts DO
  1144.       BEGIN
  1145.          xj := x^[j];
  1146.          yj := y^[j];
  1147.          IF (mode = 0) OR (yj = zero)  THEN
  1148.              wt := one
  1149.          ELSE
  1150.              wt := one/(SQR(sdy^[j]));
  1151.          xterm := wt;
  1152.          FOR n := uno TO nmax DO
  1153.          BEGIN
  1154.            sumx[n] := sumx[n] + xterm;
  1155.            xterm   := xterm * xj;
  1156.          END;
  1157.          yterm := wt * yj;
  1158.          FOR n := uno TO nterms  DO
  1159.          BEGIN
  1160.             sumy[n] := sumy[n] + yterm;
  1161.             yterm   := yterm * xj;
  1162.          END;
  1163.       END;
  1164.  
  1165. { determine the regression coefficients }
  1166.       FOR j := uno TO nterms DO
  1167.          FOR k := uno TO nterms DO
  1168.          BEGIN
  1169.            n         := j + k - uno;
  1170.            arry[j,k] := sumx[n];
  1171.          END;
  1172.        delta := determ(arry,nterms);
  1173.       IF delta = zero THEN
  1174.          bad_data(1)
  1175.       ELSE
  1176.       BEGIN
  1177.          FOR l := uno TO nterms DO
  1178.          BEGIN
  1179.             FOR j := uno TO nterms DO
  1180.             BEGIN
  1181.                FOR k := uno TO nterms DO
  1182.                BEGIN
  1183.                  n := j + k - uno;
  1184.                  arry[j,k] := sumx[n];
  1185.                END;
  1186.                arry[j,l] := sumy[j];
  1187.             END;
  1188.             a^[l] := determ(arry,nterms) / delta
  1189.          END;
  1190.  
  1191. { calculate r and standard error of estimate }
  1192.          srs    := zero;
  1193.          sum_y  := zero;
  1194.          sum_y2 := zero;
  1195.          FOR j := uno TO npts DO
  1196.          BEGIN
  1197.            yj      := y^[j];
  1198.            xj      := x^[j];
  1199.            ycalc   := zero;
  1200.            x_power := one;
  1201.            FOR k := uno TO nterms DO
  1202.            BEGIN
  1203.                ycalc   := ycalc + a^[k] * x_power;
  1204.                x_power := x_power * xj;
  1205.            END;
  1206.            residual := ycalc - yj;
  1207.            srs      := srs + SQR(residual);
  1208.            sum_y    := sum_y + yj;
  1209.            sum_y2   := sum_y2 + (yj * yj);
  1210.          END;
  1211.  
  1212.          sigmay2 := (sum_y2 - SQR(sum_y)/npts)/(npts-1);
  1213.          IF sigmay2 <= zero THEN
  1214.             bad_data(2)
  1215.          ELSE
  1216.          BEGIN
  1217.               se2    := srs/(npts - 2);
  1218.               dummy  := se2/sigmay2;
  1219.               IF  dummy <= one THEN
  1220.               BEGIN
  1221.                   r := SQRT(one - dummy);
  1222.                   IF (nterms = 2) AND (a^[nterms] < zero) THEN
  1223.                      r := -r;
  1224.                   se := SQRT(se2);
  1225.               END
  1226.               ELSE
  1227.                   bad_data(3);
  1228.          END;
  1229.       END;
  1230. END;
  1231.  
  1232. PROCEDURE smooth121(n:WORD; VAR y: single_array_pointer);
  1233. { Author : Norton Associates
  1234.   Purpose: Smooth data via binomial 121
  1235.   Version: 1.0
  1236.   Date   : 5 May 1990 }
  1237.  
  1238. VAR
  1239.    nmax,j : WORD;
  1240.    y_old,
  1241.    new_y  : SINGLE;
  1242.  
  1243. BEGIN
  1244.  
  1245. { check to see if there is enough data }
  1246.      IF n < 3 THEN
  1247.      BEGIN
  1248.          WRITELN('smooth121> no sample data');
  1249.          EXIT;
  1250.      END;
  1251.      nmax  := n - 1;
  1252.      y_old :=  y^[1];
  1253.  
  1254. { smooth data }
  1255.      FOR j := 1 TO nmax DO
  1256.      BEGIN
  1257.           new_y :=  (y_old + c2 * y^[j] + y^[j+1]) * i_4;
  1258.           y_old := y^[j];
  1259.           y^[j] := new_y;
  1260.      END;
  1261. { edge effect }
  1262.      y^[n] := (y_old + 3.0 * y^[n])* i_4;
  1263. END;
  1264.  
  1265. PROCEDURE smooth14641(n:WORD; VAR y: single_array_pointer);
  1266. { Author : Norton Associates
  1267.   Purpose: Smooth data via binomial 14641
  1268.   Version: 1.0
  1269.   Date   : 5 May 1990 }
  1270.  
  1271. VAR
  1272.    nmax,j  : WORD;
  1273.    y_2,y_1,
  1274.    new_y   : SINGLE;
  1275.  
  1276. BEGIN
  1277.  
  1278. { see if there is enough data }
  1279.      IF n < 5 THEN
  1280.      BEGIN
  1281.          WRITELN('smooth14641> no sample data');
  1282.          EXIT;
  1283.      END;
  1284.      nmax  := n - 2;
  1285.      y_2   := y^[1];
  1286.      y_1   := y^[1];
  1287.      FOR j := 1 TO nmax DO
  1288.      BEGIN
  1289.           new_y := (y_2 + c4*y_1 + c6*y^[j] + c4*y^[j+1] + y^[j+2])* i_16;
  1290.           y_2   := y_1;
  1291.           y_1   := y^[j];
  1292.           y^[j] := new_y;
  1293.      END;
  1294.      new_y := y^[n-1];
  1295.  
  1296. { edge effects }
  1297.      y^[n-1] := (y_2 + c4*y_1*  c6*y^[n-1] + c5*y^[n]) * i_16;
  1298.      y^[n]   := (y_1 + c4*new_y + c11*y^[n]) * i_16;
  1299. END;
  1300.  
  1301. PROCEDURE smoothcurve(n:WORD; VAR y: single_array_pointer);
  1302. { Author : Norton Associates
  1303.   Purpose: Smooth curvelinear data
  1304.   Version: 1.0
  1305.   Date   : 5 May 1990 }
  1306.  
  1307. VAR
  1308.    nmax,j    : WORD;
  1309.    y_2,y_1,
  1310.    new_y     : SINGLE;
  1311.  
  1312. BEGIN
  1313.  
  1314. { check to see if there is enough data }
  1315.      IF n < 6 THEN
  1316.      BEGIN
  1317.          WRITELN('smoothcurve> no sample data');
  1318.          EXIT;
  1319.      END;
  1320.      nmax := n - 2;
  1321.  
  1322.      y_2 := y^[1];
  1323.      y_1 := y^[1];
  1324.      FOR j := 1 TO nmax DO
  1325.      BEGIN
  1326.           new_y := ( c3*(y_2 + y^[j+2]) + c12*(y_1 + y^[j+1]) + c17*y^[j])* i_35;
  1327.           y_2   := y_1;
  1328.           y_1   := y^[j];
  1329.           y^[j] := new_y;
  1330.      END;
  1331.      new_y := y^[n-1];
  1332.  
  1333. { edge effects }
  1334.      y^[n-1] := (y_2 + c4*y_1 * c6*y^[n-1] + c5*y^[n]) * i_16;
  1335.      y^[n]   := (y_1 + c4*new_y + c11*y^[n]) * i_16;
  1336. END;
  1337.  
  1338. PROCEDURE corcoef(n: WORD; x,y:single_array_pointer; VAR r:SINGLE);
  1339. { Author : Norton Associates
  1340.   Purpose: Calculate correlation coefficient based upon
  1341.            product moment algrothim
  1342.   Version: 1.0
  1343.   Date   : 5 May 1990 }
  1344.  
  1345. VAR
  1346.    sumx,sumy,sumxy,sumy2,sumx2 : DOUBLE;
  1347.    j : WORD;
  1348. BEGIN
  1349.  
  1350. { check to see if there is enough data }
  1351.     IF n < dos THEN
  1352.      BEGIN
  1353.           WRITELN('corcoef> no sample data');
  1354.           EXIT;
  1355.      END;
  1356. { zero out some data }
  1357.      sumx  := 0.0;
  1358.      sumy  := 0.0;
  1359.      sumxy := 0.0;
  1360.      sumy2 := 0.0;
  1361.      sumx2 := 0.0;
  1362.  
  1363. { calculate sum of squares }
  1364.      FOR j := 1 TO n DO
  1365.      BEGIN
  1366.           sumx  := sumx + x^[j];
  1367.           sumy  := sumy + y^[j];
  1368.           sumxy := sumxy + x^[j]*y^[j];
  1369.           sumx2 := sumx2 + SQR(x^[j]);
  1370.           sumy2 := sumy2 + SQR(y^[j]);
  1371.      END;
  1372.  
  1373. {calculate correlation coefficient }
  1374.      r := ( (n*sumxy) - (sumx*sumy) )/
  1375.           ( ( SQRT((n*sumx2 - SQR(sumx))) * (SQRT(n*sumy2 - SQR(sumy))) ));
  1376. END;
  1377.  
  1378. PROCEDURE mulreg(n : WORD;
  1379.              y,x,z : single_array_pointer;
  1380.              VAR a : single_array_pointer;
  1381.              VAR r : SINGLE;
  1382.              VAR se: SINGLE);
  1383. { Author : Norton Associates
  1384.   Purpose: Determine via least squares the regresion
  1385.            and correlation coefficients for multiple regression
  1386.   Version: 1.0
  1387.   Date   : 5 May 1990 }
  1388.  
  1389. VAR
  1390.   smx,smy,smz,smz2,smx2,smxz,smxy,smyz,
  1391.   delta,srs,sum_y,sum_y2,
  1392.   ycalc,x_power,residual : EXTENDED;
  1393.   sumy                   : ARRAY[1..maxorder] OF EXTENDED;
  1394.   arry                   : arry_type;
  1395.   j,nmax,k,l             : WORD;
  1396.   xj,yj,zj               : SINGLE;
  1397.   sigmay2,se2,dummy      : EXTENDED;
  1398.  
  1399. PROCEDURE bad_mul_data(nn : WORD);
  1400. VAR
  1401.    j : WORD;
  1402. BEGIN
  1403.      r      := zero;
  1404.      se     := zero;
  1405.      FOR j := 1 TO 3 DO
  1406.          a^[j] := zero;
  1407.      WRITELN('mulreg> Determinate or correlation = 0.0   ',nn);
  1408. END;
  1409.  
  1410. PROCEDURE coeff(VAR coef : SINGLE; kk : WORD ; tarry : arry_type);
  1411. VAR
  1412.    i : WORD;
  1413. BEGIN
  1414.      FOR i := 1 TO 3 DO
  1415.      BEGIN
  1416.           tarry[i,kk] := sumy[i];
  1417.           IF kk > 1 THEN
  1418.              tarry[i,kk - 1] := arry[i,kk - 1];
  1419.      END;
  1420.      coef := determ(tarry,3) / delta;
  1421. END;
  1422.  
  1423. BEGIN
  1424.  
  1425. { check to see if there is enough data }
  1426.       IF n <= 3 THEN
  1427.       BEGIN
  1428.          WRITELN('mulreg> sample data less than 4');
  1429.          EXIT;
  1430.       END;
  1431. { zero data out }
  1432.       r := zero;
  1433.       FOR j := 1 TO 3 DO
  1434.           a^[j] := zero;
  1435.       smy   := zero;
  1436.       smz   := zero;
  1437.       smx   := zero;
  1438.       smx2  := zero;
  1439.       smz2  := zero;
  1440.       smxz  := zero;
  1441.       smyz  := zero;
  1442.       smxy  := zero;
  1443.  
  1444. { calculate sum of squares for data }
  1445.       FOR j := 1 TO n DO
  1446.       BEGIN
  1447.               yj   := y^[j];
  1448.               xj   := x^[j];
  1449.               zj   := z^[j];
  1450.               smy  := smy + yj;
  1451.               smx  := smx + xj;
  1452.               smz  := smz + zj;
  1453.               smx2 := smx2 + SQR(xj);
  1454.               smz2 := smz2 + SQR(zj);
  1455.               smxz := smxz + (xj*zj);
  1456.               smxy := smxy + (xj*yj);
  1457.               smyz := smyz + (yj*zj);
  1458.       END;
  1459.  
  1460. { set up linear simultaneous equations}
  1461.       sumy[1]   := smy;
  1462.       sumy[2]   := smxy;
  1463.       sumy[3]   := smyz;
  1464.       arry[1,1] := n;
  1465.       arry[2,1] := smx;
  1466.       arry[3,1] := smz;
  1467.       arry[1,2] := smx;
  1468.       arry[2,2] := smx2;
  1469.       arry[3,2] := smxz;
  1470.       arry[1,3] := smz;
  1471.       arry[2,3] := smxz;
  1472.       arry[3,3] := smz2;
  1473.  
  1474.       delta := determ(arry,3);
  1475.       IF delta = zero THEN
  1476.          bad_mul_data(1)
  1477.       ELSE
  1478.       BEGIN
  1479.          coeff(a^[1],1,arry);
  1480.          coeff(a^[2],2,arry);
  1481.          coeff(a^[3],3,arry);
  1482.       END;
  1483.  
  1484. { calculate r and standard error of estimate }
  1485.       srs    := zero;
  1486.       sum_y  := zero;
  1487.       sum_y2 := zero;
  1488.       FOR j := 1 TO n DO
  1489.       BEGIN
  1490.            yj       := y^[j];
  1491.            xj       := x^[j];
  1492.            zj       := z^[j];
  1493.            ycalc    := a^[1]  + a^[2]*xj + a^[3]*zj;
  1494.            residual := ycalc - yj;
  1495.            srs      := srs + SQR(residual);
  1496.            sum_y    := sum_y + yj;
  1497.            sum_y2   := sum_y2 + SQR(yj);
  1498.       END;
  1499.  
  1500.       sigmay2 := (sum_y2 - SQR(sum_y)/n)/(n - 1);
  1501.       IF sigmay2 <= zero THEN
  1502.          bad_mul_data(2)
  1503.       ELSE
  1504.       BEGIN
  1505.            se2    := srs/(n - 3);
  1506.            dummy  := se2/sigmay2;
  1507.            IF  dummy <= one THEN
  1508.            BEGIN
  1509.               r  := SQRT(one - dummy);
  1510.               se := SQRT(se2);
  1511.            END
  1512.            ELSE
  1513.               bad_mul_data(3);
  1514.      END;
  1515. END;
  1516.  
  1517. PROCEDURE autocor(n:WORD; x:single_array_pointer; lag :WORD; VAR auto : single_array_pointer);
  1518. { Author : Norton Associates
  1519.   Purpose: Determine correlation coefficients for auto lagging of the data
  1520.   Version: 1.0
  1521.   Date   : 5 May 1990 }
  1522.  
  1523. VAR
  1524.    k,j,
  1525.    num  : WORD;
  1526.    y    : single_array_pointer;
  1527.    r    : SINGLE;
  1528.  
  1529. BEGIN
  1530.  
  1531. { check to see if there is enough data }
  1532.     IF n < dos THEN
  1533.     BEGIN
  1534.          WRITELN('autocor> no sample data');
  1535.          EXIT;
  1536.     END;
  1537.  
  1538. { lag value can not be less than the number of data points }
  1539.     IF lag > n THEN
  1540.     BEGIN
  1541.          WRITELN('autocor> lag can not exceed number in sample data');
  1542.          EXIT;
  1543.     END;
  1544.     create_single_array(n,y);
  1545.     FOR j := 1 TO lag DO
  1546.     BEGIN
  1547.          num := n - j + 1;
  1548.          MOVE(x^[j],y^[1],(num*SIZEOF(x^[1])));
  1549.          corcoef(num,x,y,r);
  1550.          auto^[j] := r;
  1551.     END;
  1552.     delete_single_array(n,y);
  1553.  
  1554. END;
  1555.  
  1556. FUNCTION movavg(n: WORD; a : single_array_pointer; ma:WORD; k : WORD) : SINGLE;
  1557. { Author : Norton Associates
  1558.   Purpose: Determine the moving average of data for a selected data set
  1559.   Version: 1.0
  1560.   Date   : 5 May 1990 }
  1561.  
  1562. VAR
  1563.    dum : SINGLE;
  1564.    j   : WORD;
  1565.  
  1566. BEGIN
  1567. { check to see if index is out of range }
  1568.      IF (k > n) OR ( k < 1) THEN
  1569.      BEGIN
  1570.           WRITELN('movavg> sorry index is out of range');
  1571.           EXIT;
  1572.      END;
  1573. { check to see if the index is less than the lag number }
  1574.      IF k < ma THEN
  1575.      BEGIN
  1576.           WRITELN('movavg> number of points less than moving average: movavg = 0');
  1577.           movavg := 0.0;
  1578.      END
  1579.      ELSE
  1580.      BEGIN
  1581.          dum := 0.0;
  1582.          FOR j := (k - ma + uno) TO k DO
  1583.              dum := dum + a^[j];
  1584.          movavg := dum  /  ma ;
  1585.      END;
  1586. END;
  1587. BEGIN
  1588.  
  1589. END.
  1590.